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We study the nonlinear rheology of a glass-forming binary mixture under the reversal of shear flow using 
molecular dynamics simulations and a schematic model of the mode-coupling theory of the glass transition 
(MCT). Memory effects lead to a history-dependent response, as exemplified by the vanishing of a stress- 
overshoot phenomenon in the stress-strain curves of the sheared liquid, and a change in the apparent clastic 
coefficients around states with zero stress. We investigate the various retarded contributions to the stress 
response at a given time schematically within MCT. The connection of this macroscopic response to single- 
particle motion is demonstrated using molecular-dynamics simulation. 



I. INTRODUCTION 

Dense liquids exhibit slow relaxation processes in ther- 
mal equilibrium, characterized by a time scale r that di- 
verges as the liquid is driven to structural arrest into 
an amorphous solid^. This opens a window rich in non- 
linear phenomena under flow, both amenable to experi- 
ment and theoretical treatment. It is of interest, because 
the externally imposed flow strongly disturbs the equilib- 
rium relaxation processes. More precisely, given a typical 
single-particle relaxation time scale r , we consider shear 
flow with deformation rate 7 such that Pe$ = 7T0 <C 1, 
but Pe = 7T ^> 1. This nonlinear rheology has in partic- 
ular been studied in colloidal dispersions and soft-matter 
systems, where the mesoscopic size of the constituent par- 
ticles implies relaxation times of 0(s), easily probed by 
moderate flow in experiment. In the flowing steady state, 
shear thinning is the most prominent phenomenon;^ the 
apparent viscosity of the sheared fluid decreases rapidly 
with increasing shear rate 7. Quite generally, the flow- 
ing steady state is characterized by a nonlinear relation 
between the macroscopic stress a and the applied shear 
rate 7. 

Going beyond the steady state, the behavior under 
start-up flow has received particular attention. As a 
steady shear flow of rate 7 is suddenly switched on at 
t = 0, stresses a increase from zero (assuming the exper- 
iment was started in an equilibrated and hence stress- free 
configuration) to their steady-state value (7(7). The re- 
sulting stress-strain curve (7(7), measured as a function 
of applied deformation 7 = <yb for t > 0, exhibits a pe- 
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culiar intermediate maximum, termed stress overshoot. 
This is typically found in amorphous metals^, or (for dif- 
ferent reasons) polymer melts^. The maximum of the 
(7(7) curve has been linked to material stability: it can 
be taken as a yield stress that separates reversible (linear 
and nonlinear) elastic deformation at small strains from 
irreversible plastic flow induced by larger deformations. 
More recently, stress overshoots have also been discussed 
in soft colloidal systems, underlining a generic mechanism 
applicable in a wide variety of system classes^—. 

As a transient phenomenon, the stress overshoot is sus- 
ceptible to the history of sample preparation. Molecular 
dynamics simulations have addressed this focusing on the 
shear-molten amorphous solid, where a systematic de- 
pendence was found on the waiting time t age between 
sample preparation and flow start-up^, as long as the 
start-up configuration is not fully equilibrated, t agc <C r. 

In the engineering literature it is well established that 
the stress-strain curve depends on the deformation his- 
tory of the sample. One particular example was found 
empirically by Johann Bauschinger in the 1880's during 
experiments on various steels^: these materials typically 
are softer under compressional load after having been 
subject to tension. This is, in the strict sense, called 
the Bauschinger effect. As an abstraction, take the case 
where pre-shear is applied in one particular direction, 
and consider the modified stress-strain curve measured 
if shear is then applied in the opposite direction. Typi- 
cally, the regime of elastic deformation shrinks, so that 
plastic rearrangements dominate much earlier than in an 
initially unstrained sample. 

Reversing the shear rate, which we will call 
'Bauschinger effect' also, in the following, has 
been discussed^ 2 - in the framework of the shcar- 
transformation-zone theory (STZ), where one attributes 
internal degrees of freedom to the STZs, which carry the 
memory of the deformation history. Recent simulation 
studies of the Bauschinger effect^ 3 - in a system quenched 
into its low-temperature glassy state have been analyzed 
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in terms of anisotropic elastic constants that arise due to 
the deformation history. 

It is a matter of current debate whether and how 
the rheological properties of amorphous materials deep 
in the glassy state (here referred to as essentially the 
T — > case), and at much higher temperaturs in the melt, 
are related. In athermal systems, concepts like STZs^ 
and avalanches of yielding events^ are widely recognized 
as useful to rationalize macroscopic behavior. But the 
range of applicability of these approaches is not well un- 
derstood when dealing with liquids close to their glass 
transition, where thermal fluctuations dominate yield- 
ing mechanisms that may be qualitatively different from, 
say, athermal jamming^. The rheology of dense liquids 
is the realm of a recent extension of the mode-coupling 
theory of the glass transition (MCT) to sheared colloidal 
systems^—. Intriguingly, some of the macroscopic phe- 
nomena are quite alike in both quasi-athcrmal systems 
and liquids driven by thermal noise. 

Here we investigate, both by mode-coupling theory and 
computer simulation, the nonlinear rheology under flow 
reversal in the initially equilibrated liquid state, without 
referring to T — > properties. The initial stress over- 
shoot is shown to vanish if the flow is reversed in the 
steady state. Hence, the Bauschinger effect does not re- 
quire an explanation associated with properties of the 
quenched glass or aging phenomena. We study its con- 
nection to the time-dependent memory of past flow that 
is typical for dense liquids, by investigating the effect 
of flow reversal applied at various times during the ini- 
tial evolution from equilibrium to steady state. A clear 
separation between reversible elastic rearrangements and 
irreversible plastic flow emerges. We connect the macro- 
scopic stress response to microscopic particle motion, 
quantified through the transient mean-squared displace- 
ments. 

The paper is organized as follows: Section Ull contains 
the exposition of the theoretical approach, and section HTIl 
gives an overview of the simulation. In Sec. [TV] we discuss 
the comparison between theory and simulation, and the 
insights thereby gained. Section [V] concludes with an 
outlook. 



II. THEORY 

Our theoretical approach is based on the mode- 
coupling theory of the glass transition (MCT)iS, as ex- 
tended to deal with nonlinear rheology in the integration- 
through transients (ITT) formalism!^. ITT-MCT is a 
microscopic (albeit approximate) theory describing the 
dynamical correlation functions of the wave-vector de- 
pendent density fluctuations in glass-forming liquids. 
Suitably generalized Green-Kubo relations provide the 
nonequilibrium transport coefficients in the nonlinear- 
response regime; within MCT, they are expressed as in- 
tegrals over the density correlators. The coupling coef- 
ficients of the theory arc fully determined by knowing 



the equilibrium static structure of the liquid. The the- 
ory then predicts an ideal glass transition to occur as 
a bifurcation transition: smooth changes in control pa- 
rameters such as density or temperature bring about a 
discontinuous change in the long-time behavior of the 
density correlation functions. These decay to zero in the 
equilibrium liquid, but attain a finite long-time limit, the 
glass form factor, in the ideal glass. The theory compares 
semi-quantitativcly to the linear rheology measured in 
colloidal hard sphere dispersions 2 ^— . 

Since the full microscopic theory is rather cumber- 
some to solve, various simplifications have been devised. 
Schematic models try to capture mathematically the 
essence of the bifurcation transition by reducing the num- 
ber of correlation functions to a few or only one, at the 
cost of neglecting spatial information. The important 
correlations contained in the schematic models are tem- 
poral ones. Memory effects are contained in general- 
ized friction kernels depending on (up to) three times. 
Causality is automatically enforced, and the flow history 
of the sample is contained in the friction-kernels via time- 
or strain-dependent clastic coefficients and structural re- 
laxation functions. The latter describe plastic deforma- 
tion, and their nonlinear equations of motion are the 
heart of the theory. The non-Markovian equations of the 
schematic models contain a small number of fitting pa- 
rameters that mimic the relation between the equilibrium 
structure (encoding the interaction potential) and the 
theory's coupling coefficients. Smooth and regular vari- 
ations of the parameters lead to dramatic elasto-plastic 
and viscoelastic changes in the structural relaxation and 
non-linear rheology. The interpretation of the effects in 
terms of spatially dependent particle rearrangements can 
not be achieved by the schematic model, and thus it re- 
mains unresolved to which degree 'force-chains', 'hetero- 
geneities', or specific 'plastic rearrangement centers' are 
described in a (presumably) space-averaged manner. 



A. Schematic MCT Model 

A schematic MCT model for the time-dependent rhe- 
ology of colloidal suspensions was introduced in Ref. [23[ ; 
it has been used to study large-amplitude oscillatory 
shear—, mixed flows2&, and double-step strain 2 ^, and 
generalizes a model used to analyze stationary flow 
curve s 20 ' 21 . In the model, the nonlinear generalization 
of the Green-Kubo relation for the shear-stress tensor el- 
ement a = a xy reads, imposing spatially homogeneous 
but otherwise arbitrary simple shear flow with rate 7 in 
the ^-direction, with gradient in the y-direction, 

a{t)= J dt'j(t')G(t,t',m. (1) 

—00 

The generalized shear modulus G(t, t', [7]) depends on 
the full flow history encoded in the time-dependent shear 
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rate j(t). Outside the steady state, it also depends on 
the two times corresponding to the underlying fluctua- 
tions separately. In the spirit of MCT, we approximate 
the shear modulus in terms of the transient density cor- 
relation function $(f, t'), 

G(t,t')=v a (t,t')<S> 2 (t,t'). (2) 

Here, an elastic coefficient v a appears that was set to a 
constant in the original formulation of the model^. It 
measures the strength of stress fluctuations caused by 
pair-density fluctuations. The microscopic theory im- 
plies a relation between &(t,t') and G(t,t') that involves 
an integration over wave vectors, including nontrivial 
weights that depend on time since density fluctuations 
are advected by shear. It is therefore plausible to extend 
the original schematic model by allowing for a prefactor 
v a (t,t') that depends on time through the accumulated 
strain j(t, t') = L j(s) ds, 

*M-<f(^)X-(^)1- 

(3) 

Symmetry dictates that the direction of the strain does 
not matter, leading to an even function v a — v a {^(t, t') 2 ). 
Equation © generalizes the time-dependent coupling 
first introduced for constant shear rates in Ref. @ , where 
it was justified by a comparison with fully microscopic 
ITT-MCT calculations in two dimensions. The func- 
tional form is chosen simple enough for rapid numerical 
implementation and allows for a quantitative discussion 
of stress overshoots in model colloidal suspensions^* 2 !. 
The parameter v* sets the scale of the generalized shear 
modulus. We obtain it by matching it to the linear- 
response low-frequency plateau modulus G M of the sim- 
ulation. The shear elastic constant G^ itself often is 
called 'shear modulus', and can, e.g., be measured in 
the low-frequency linear elastic modulus of the glass, 
G'(co 0) = Gqo, or as one of the two Lame-coefficients 
of an isotropic solid^. At a critical strain scale 7*, 
the time-dependent modulus becomes negative and the 
stress-strain curve exhibits an overshoot with its peak at 
the critical strain. 7** describes the decay of the over- 
shoot. To reduce the number of free parameters, we keep 
it linearly proportional^ to 7* . 

The time-dependent elastic coefficient given in Eq. ([3]) 
is carefully chosen to become negative in a small strain 
window relevant for the decay of the correlation func- 
tions. Doing so, it allows to circumvent an obvious limita- 
tion of the schematic model: setting v„ = const., Eq. @ 
has a definite sign. In startup shear, cf. Eq. ([1]), this 
implies that a(t) varies monotonically with accumulated 
strain. The original schematic model hence does not de- 
scribe stress overshoot phenomena. Motivated by the full 
microscopic expression of ITT-MCT—^ and the obser- 
vation that the wave-vector dependent MCT, even with 
further isotropy assumptions, produces a stress overshoot 
through a small negative dip in G(t,t'), cf. Ref. Q, 



Eq. Q is a simple way to incorporate the phenomenol- 
ogy. In the following, we take it as one established way 
to match the startup stress-strain curves (see below for a 
comparison to simulation), and discuss the consequences 
of the history-integral formulation of Eq. (JTJ) under flow 
reversal. 

The schematic density correlator is the solution of a 
Mori-Zwanzig-likc memory equation 

<b{t,t') + T ^>(t,t') + jdt"m(t,t",t')^(t",t')j = 0, 

(4) 

with initial condition t') = 1 - T{t — t') + . . .. Here 
r denotes the initial decay rate of the density correlator, 
which corresponds to a reciprocal microscopic relaxation 
time. We take T to be independent of the shear rate. 

To close this equation of motion, MCT approximates 
the memory kernel m through a nonlinear polynomial 
of the correlation functions. In contrast to a Marko- 
vian approximation, this postulates that the kernel re- 
laxes on the same time scale as the correlator, and that 
both need to be obtained self-consistently. Conceptually, 
this is motivated by the assumption that stress fluctua- 
tions captured in m arise from slow structural rearrange- 
ments, which themselves are captured in $. The simplest 
schematic model that recovers the asymptotic behavior 
found for typical quiescent glass-forming liquids is the 
F12 models, whose extension to time-dependent shear 
reads 

m(t, t", t') = h(t, t") [«!$(£, t") + v 2 $ 2 (t, t")] . (5) 

The parameters v\ and v 2 are the coupling coefficients 
driving the system through the glass transition. There is 
a line of bifurcation points (vf, v 2 ), and we pick v 2 = v% 
and vi = vl + e/(V2 - 1) with v{ = 2{s/2 - 1) as in 
previous work. The parameter e controls the distance to 
the glass transition: e < corresponds to liquid states, 
where the correlation function decays to zero in the qui- 
escent system, on a time scale r that diverges as e — > 0. 
States with e > are idealized glass states, where the qui- 
escent correlation function attains a positive long-time 
limit, lirtii-^oo^^t' <f>(t,i') = / > 0. This value is called 
the non-ergodicity parameter or glass form factor. The 
clastic constant in the model then follows as Goo = v%f 2 . 

One of the predictions of ITT-MCT is that (steady) 
shear melts the glass; the correlation functions even in 
the glass show shear-induced decay due to a loss of 
memor y 16 ' 29 . In the schematic model, this is modeled 
by an ad-hoc strain-reduction function, 

h(t,t') = ^. (6) 

1 + ( 7 (M')/7c) 2 

Symmetry again dictates that h is an even function of 
the accumulated strain. The parameter 7 C is a critical 
strain and describes how fast the memory of the glass de- 
cays due to its deformation. In principle, the microscopic 
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Figure 1. Linear plot of the stress-strain relation (7(7) at fixed 
density and shear rate for various flow histories: startup flow 
with constant shear rate (red), flow reversal in the steady 
state (s, purple), from the elastic regime (el, blue), and from 
the point of the stress overshoot (max, green). The linear in- 
crease with the quiescent elastic modulus Goo in startup, and 
with G e ff after shear-reversal are indicated by dotted lines. 
Characteristic stress-strain points, whose transient modulus 
is discussed in the text, are marked by symbols: The red 
square gives the stationary stress. It corresponds to an inte- 
gration along the line from A to B in Fig. [7] Triangles mark 
stresses after flow reversal in the linear regime and correspond 
to an integration along the line C to E'. Discs mark stresses 
after flow reversal in the overshoot regime and are obtained 
via integrations along the line C to E in Fig. 

form of the MCT equation suggests further combinations 
of the three times to appear in Eq. (|5|), as exploited for 
example in Refs. [23[ and (2(|. We neglect these terms 
here, to make contact with earlier schematic-model anal- 
yses of the nonlinear steady-state rheology of colloidal 
suspensions^— . In principle, the strain-reduction func- 
tion h could also contain negative regions, as arises in the 
stress-density coupling v a . This could cause the transient 
correlator $ to become negative for intermediate times as 
seen in computer simulations and microscopic ITT-MCT 
calculations in two dimensions^. For simplicity, this is 
neglected. 

B. Bauschinger effect 

The flow history we consider in this work is given by 
two intervals of constant shear rate: 

!0 for t < 0, 
+7 for0<i<i w , (7) 
—7 for t > t w . 

We mainly discuss a fixed value for 7, and different wait- 
ing times t w over which the "pre-shear" in the positive 
direction is applied. These times correspond to an accu- 
mulated strain 7^ = jt w . The amount of this pre-strain 
is decisive and determines the stress-strain curve after 



switching. Figure [TJ presents typical results obtained in 
the schematic model as unmodified stress-versus-strain 
curves; it serves to introduce the flow history and char- 
acteristic parameters. After startup of shear at t = 0, the 
system first undergoes an elastic transient, characterized 
by a stress-strain curve 17(7) that is approximately linear, 
cr w Goo7, given by the quiescent elastic shear modulus. 
At the end of this regime, the cr-versus-7 curves becomes 
sublincar, until the position of the overshoot is reached 
at some strain 7 max ; typically of the order of 10%. At 
large strains, one reaches the steady state characterized 
by a constant (7(7 — > 00) = a s independent of 7. For 
the flow reversal, we consider switching the shear rate 
at three typical points in either regime: 7^ inside the 
elastic transient, 7™ ax , and 7^ in the steady state. The 
7k, values are marked in Fig. [TJ A value 7^, = 0.4 in 
the schematic model turns out to be sufficient so that 
no further "f w dependence is observed for larger values 
besides a trivial shift. To ease the interpretation of the 
curves at different switching times, we choose 7^ such 
that the corresponding stress is equal to the steady-state 
one, (7(7^) = ct(t'5j) ■ Generally, after reversing the flow 
direction, the stress first decreases back to zero, defining 
a strain value 70. The stress-strain curves in later figures 
will be compared in a | a\ -versus- 1 7 — 70I representation. 
While it does hide the "hysteresis" loop which can be rec- 
ognized in Fig. [TJ and which is familiar in the engineering 
literature, it simplifies the quantitative analysis. 

All calculations in the schematic model are done with 
parameters that have been found typical in applications 
of predecessors of the model to large amplitude oscilla- 
tory shear—, and especially to startup flow^i. The val- 
ues are: T = 100, j c = 0.75, 7* = 0.1, and 7** = 0.133; 
see also Fig. [2] The considered state in Fig. [TJ which 
will be analyzed in more detail in Figs. [5] and [5J is a 
glass very close to the glass transition in order to achieve 
a large separation of the structural dynamics from the 
short-time one; e = 10~ 4 . The shear rate 7 = 5 x 10~ 3 , 
for which the Bauschinger effect is studied, corresponds 
to a rather small Peclet number, which in the model is 
given by Pe = 7/T, and thus takes the value 5 x 10~ 5 . 
These two values provide access to the asymptotic regime 
of ITT-MCT for e ^ and Pe < 1. Slightly different e 
and 7 values are used in the comparisons to the simula- 
tions in Sec. IIVI 

FigurcUJcontains a wealth of information, which will be 
discussed in Sec. IIVI together with the simulation result s. 



III. SIMULATION 

We perform nonequilibrium molecular dynamics com- 
puter simulation for a glass-forming binary mixture. Par- 
ticles interact through a purely repulsive soft-sphere po- 
tential. Denoting particle species by a, /3, the truncated 
Lennard- Jones potential due to Weeks, Chandler, and 
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Andersen^ reads 



t/WCA 

V a/3 
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else 



where r r 



2 1//6 c Q ^ is a cutoff for the interaction 



range. To ensure continuity of force and conservation of 
total energy in the NVE ensemble, we apply a smooth- 
ing function, V a p(r) = V™ CA (r)S a p(r) with S a p(r) = 
(r - r c , aP ) A /[h A + (r - r c , Q/3 ) 4 ] with h = lir 2 a a ^. We 
choose units of energy such that = 1, and units of 
length such that <taa = 1- The unit of time is given 
by yj TOAC AA /eAA where to a = = 1 are the masses 
of the particles. For the smaller particles, obb = 5/6, 
and the mixture is additive, gab = (o"aa + cbb)/2. This 
system, albeit with differing masses, has been previously 
studied in the quiescent state and identified as a model 
glass former in its equimolar composition 3 ^. 

The system is coupled to a dissipative particle dynam- 
ics (DPD) thermostat: the equations of motion read^i 



Pa 



f a = F^m) + F»({r,p\) + ?*({?}) 



(8b) 



where r and p are positions and momenta of the particles, 
and F l denotes conservative (i = C), dissipative (i = 
D), and random (i = R) forces. With the interparticle 
separation = |r, — ?j\ between particles i and j of 
species a and fi, 



F2 = -J2VV a p(r 



3& 



F Q R = ^V^fcBTC^(^)AA y 



(9) 
(10) 

(11) 



Here, Vij = Pi/mi ~pj/rrij is the relative velocity be- 
tween the two particles, and rij = (fi — Tj)/rij their unit 
separation. The appearance of the relative velocity in 
the dissipative force is crucial to obtain local momentum 
conservation and Galilean invarianco^. w(r) is a cut- 
off function, set to w = 1 for r < t" c ,dpd/2 and w = 
elsewhere. We have chosen r C DPD = l-7r CjCtQ in our sim- 
ulations. C is a parameter controlling the strength of 
friction forces; we set £ = 10. This parameter is not cru- 
cial for the results to be discussed^. The j\f are Gaussian 
normal random variables. 

The DPD equations of motion, Eq. ([8]) , are integrated 
with a generalized velocity Verlet algorithm 3 ^. This al- 
gorithm integrates the reversible Hamiltonian part of the 
equations with the velocity Verlet scheme and then par- 
tially re-equilibrates the two-particle momenta, ensuring 
a Boltzmann distribution in equilibrium. A time step of 
St = 5 x 10~ 4 (in the specified units of time) is employed. 
Two different neighbor lists are implemented, one for the 



force calculations, and one for the DPD thermostat. The 
simulation consists of TV = 2N\ = 1300 particles in a 
three-dimensional box with volume V = L? and linear 
dimension L = IOoaa, corresponding to a number den- 
sity p = 1.3. At this density, no signs of crystallization 
or phase separation were observed in the studied tem- 
perature range from T = 5 to T = 0.4. A glass transi- 
tion point is estimated as T c w 0.347 according to mode- 
coupling theory. In the following we focus on the equili- 
brated fluid at T = 0.4. Initial equilibration proceeded 
by using St = 10~ 3 , assigning new velocities every 50 in- 
tegration time steps. Equilibration was checked by the 
decay of the incoherent intermediate scattering function 
at a wave number corresponding to a typical interparti- 
cle separation q = 27t/(7aa: runs were long enough to ob- 
serve the decay of the correlation function to zero at long 
times. A set of 200 independently equilibrated configu- 
rations served as initial configurations for the production 
runs employing the DPD thermostat. 

Shear is applied in the cc-direction with a gradient in 
y-dircction and a positive shear rate 7 initially. Pla- 
nar Couette flow is imposed by periodic Lees-Edwards 
boundary conditions 3 ^: the periodic image of a particle 
leaving or entering the simulation box in y-direction is 
displaced in x-direction according to the strain ±<yL. 

The macroscopic response of the system to imposed 
shear flow is measured through the stress tensor. In the 
specified coordinate system, the dominant contribution 
at low shear rates is its xy-componcnt, given by the Kirk- 
wood formula 3 ^, 



>xy 



1 / N 



,X ± 7 



(12) 

The first term is a kinematic contribution, while the sec- 
ond is the virial contribution incorporating the nonhy- 
drodynamic forces exerted on each particle. Angular 
brackets indicate canonical averages. In the following, 
we restrict the discussion of the Bauschingcr effect to 
one, rather high, exemplary shear rate, 7 = 5 x 10~ 3 . 



IV. RESULTS AND DISCUSSION 

A. Comparison of simulation and theory 

We begin by a discussion of the properties of the 
computer-simulated system under startup of flow and 
in its steady state. Qualitatively, they reproduce 
earlier results^. Here, they serve to determine the 
schematic model parameters prior to investigation of the 
Bauschinger effect. 

Figure [5] shows the stress-strain curves (7(7) obtained 
after switch-on of steady shear at t = 0, for various 
shear rates 7 at one particular temperature in the su- 
percooled liquid. For small 7, the system behaves as a 
linear viscoclastic fluid, and the stress quickly reaches its 
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Figure 2. Plot of the stress-strain curves a(-y) at fixed tem- 
perature for various shear rates. The dashed lines show the 
MCT results for T = 100, e = -10~ 4 , 7c = 0.75, V* = 97, and 
7*,/7* = 1.33. 7, is 10 _1 ,7x 10 _2 ,6x 10 _2 ,4.5x 10" 2 ,3.5x 
10 -2 , 1.5 x 10 -2 for decreasing shear rates 7. 



steady state value, rising monotonically from zero. As 
the shear rate is increased, an intermediate overshoot in 
the cr-versus-7 curve appears. Depending on 7, the typ- 
ical strain reached at the position of the overshoot is a 
few percent; at the highest shear rates this value is about 
10%. It is in agreement with earlier observations^, and 
one can argue that this typical strain corresponds to a 
typical size of nearest-neighbor cages that are being bro- 
ken by shear—. 

From the largc-7 limit of the stress-strain curves, one 
obtains the steady-state flow curve, <r s (7)- The values ex- 
tracted from our simulation and the schematic-model fits 
exhibit shear-thinning: The stress increases sublincarly 
with the rate. Note in Fig. [5] that the values of a are 
given in simulation units, where kT = 0.4. This confirms 
that the stresses we observe are dominated by thermal 
motion, since a s = 0(kT/a^, a ). This compares well to 
recent investigations of colloidal suspensions^ ' 21 ' 27 . The 
precise numerical value is found to differ among different 
systems, as it somewhat depends on the interaction po- 
tential and mixture parameters in experimental samples. 

The solid lines in Fig. [2] arc fits with the schematic 
MCT model including the time-dependent vertex for the 
Green-Kubo relation. These fits were optimized to repro- 
duce the shape of the overshoot at the largest shear rate 
(for which the Bauschinger effect will be discussed be- 
low). Further adjustments to improve agreement at the 
lower shear rates may well be possible, but note that the 
form of the overshoot is dictated by the empirical choice 
of v a {t,t'). 

The parameters of the schematic MCT model can 
be determined by fitting the simulation data in startup 
shear. Using a fit procedure similar to the one outlined 
in Ref. Q we arrive at the values listed in the caption of 
Fig. [21 The parameters T and e determine the short- and 
long-time relaxation scales To and r in the model. They 
are fixed by a comparison of the quiescent equilibrium 
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Figure 3. Linear plot of the stress-strain relation |ff(|7 — 7o|)| 
at fixed temperature and shear rate for various flow histories: 
starting from equilibrium (EQ, red line), after flow reversal in 
the steady state (S, purple), from the elastic regime (el, blue), 
and from the point of the stress overshoot (max, green). Re- 
versal at 7^ = 0.4 is in the steady state as seen from com- 
paring with reversal at 7^ = 75 (black data). Smooth lines 
result from the schematic model. 



correlation functions at a typical intermediate wave num- 
ber corresponding to the first peak of the static structure 
factor. Linear-response fluctuations also allow to deter- 
mine v*, corresponding to the plateau modulus and 
hence setting the units of stresses in the model. The pa- 
rameters 7c, 7*, and 7** are then tuned to the startup 
rhcology data, taking the shear rate 7 nominally from 
the simulation. 

Once the parameters given in the caption of Fig. [2] 
have been fixed, the predictions of the schematic model 
regarding the modification of the stress-strain curve un- 
der flow reversal are parameter-free. 

We now turn to the discussion of the Bauschinger ef- 
fect. Figure [3] shows the main result of the computer 
simulation, viz. the stress-strain curves in the |<r|-versus- 
|7 — 7o I representation, i.e., shifted and inverted such that 
they all start from a stress-free state and display the be- 
havior as a function of the additional strain imposed on 
this configuration. Curves for various j w are shown, as 
discussed above. They all agree at large strains by con- 
struction, confirming that the steady-state stress does 
not depend on the shear history. Comparing first the 
startup curve, j w = 0, with the one for reversal in the 
steady-state, -f w =7^, the Bauschinger effect becomes 
most apparent. The two curves differ in two main as- 
pects: the linear slope at small deformations is lower, 
and the overshoot is gone in the curve corresponding to 
oppositely straining the (steadily) pre-sheared configura- 
tion. This holds for shear-reversal throughout the sta- 
tionary state, as shown by the two curves for 7^, = 0.4 
and 7^ = 75 included in Fig. [3] Both these differences 
arise because the system undergoes plastic deformations 
during pre-shcar: confining the pre-shear to small strains, 
7to = 7^, 5 so that flow reversal takes place inside the 
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Figure 4. Velocity profiles v x (y) as a function of position 
along the gradient direction y for various times after appli- 
cation of a new flow direction, (a) Starting from a quiescent 
configuration (EQ), shear in +x direction, (b) Starting from 
the resulting steady state (S) at 7J, = 0.4, shear in — x di- 
rection. Insets mark the points for which profiles are shown 
along the stress-strain curves. A total 100 independent con- 
figurations are averaged within a small time window to obtain 
the depicted graphics. 



initial elastic-deformation regime, the overshoot is main- 
tained almost unchanged, as is the effective clastic coef- 
ficient extracted from the initial rise of the stress-strain 
curve. The cross-over between elasticity-dominated and 
plasticity-dominated pre-strain occurs gradually, as the 
curve for j w = 7™ ax exemplifies. 

The solid lines in Fig. [3] represent the schematic-MCT 
model calculations. The line corresponding to 7^ = 
is a result of the fitting procedure described above. All 
the other theory curves then follow from the structure of 
the ITT-MCT equations. They describe the computer- 
simulation data extremely well, and capture both the de- 
crease of the overshoot and the decrease of the effective 
elastic coefficient. 

Before continuing the discussion, it is worth pointing 
out that none of the features we discuss appear to be 
explained by, or connected to shear banding or other in- 



homogencous flow effects. To demonstrate this, we show 
in Fig. [4] the velocity profiles obtained from the simula- 
tion at various instants in time, both under forward shear 
and the subsequently reversed shear. To improve statis- 
tics, the instantaneous velocities were averaged over all 
the configurations and over a small time window of width 
At = 0.1. Insets in the figure mark the points for which 
the velocity profiles v x (y) are shown; they all follow the 
linear behavior v x (y) = jy expected from homogeneous 
simple shear, to within statistical noise. Similar results 
have also been obtained under startup shear—. It appears 
that the small time needed to establish a linear veloc- 
ity profile under Lees-Edwards boundary conditions is 
not relevant for our discussion. In Fig. [4] only the case 
for flow reversal in the steady state is shown, but qual- 
itatively the same result is found for the velocity pro- 
files obtained after switching the shear direction at the 
smaller 7™ ax or 7™. Stress overshoot phenomena have 
been linked to shear banding^, but it appears that this 
mechanism is not relevant for our system. Indeed, our 
results are compatible with a recent study suggesting a 
density-dependent critical Pcclet number for the appear- 
ance of shear bands^S. 



B. Interpretation in a generalized Maxwell model 

To understand the history dependence of the stress- 
strain curves, let us split the integral in Eq. ([T]) into two 



parts. For t > t w , write 



o(i) = oi(t) - on(i) 



with the summands 



t w 



01 (t) =7/ dt'G(t,t'), 



(13) 



(14a) 



on(t)=7jf dt'G{t,t') =i dt'G{t-t'). (14b) 

In the last term, on(t), the condition t w < t' < t is 
fulfilled, so that the correlation and strain functions en- 
tering the integral arc dependent on one-time only and 
describe the transient between an equilibrium and a sta- 
tionary state. Hence, this contribution is trivially related 
to the startup curve, —au(t + t w ) = o s tart(i); see below 
in Sect. ITVCl for details. 

To simplify the discussion, let us consider an ad-hoc 
generalization of the Maxwell model of viscoclastic liq- 
uids to shear-thinning fluids. Keeping the form of the 
Green-Kubo equation, we approximate the generalized 
shear modulus, Eq. ([2]), by 

G(t, t') « G M (t, t') = v a (t, ^)e-('-*')l7(OI/7c _ (15) 

For the case of constant shear and constant v a , this form 
shows the same qualitative features as the shear-molten 
glass in schematic MCT— . To generalize this nonlinear 
Maxwell model to unsteady flows is neither trivial nor 



8 



(7 




Figure 5. Bauschinger effect as illustrated by a generalized 
Maxwell model (see text). Solid lines are a(t) as a function of 
normalized time, \ j\t, for the case y w — 0.1 corresponding to 
7™ ax and 7™ = 0.8 corresponding to 7^. Dashed lines display 
the contributions ai(t) (upper, red) and <rii{i) (lower, blue). 



unique, but assuming the rate of decorrelation to be set 
by the instantaneous shear rate, the above ansatz appears 
plausible, if somewhat crude^. In the present case it even 
simplifies further because |7(i')l = 7- The integrals de- 
termining a(t) and its two contributions g\ and o\\ can 
then be solved analytically when inserting the form of 
v a (t,t') given above. In order to quantitatively match 
the startup stress-strain curve of the generalized Maxwell 
model with the schematic model, we set -% = 7 c /3. 

The response for large t w can now be understood, re- 
calling that the integral appearing in a\(t) is dominated 
by t' w t w , where the integrand is still non- vanishing. 
The accumulated strain entering the correlation func- 
tions can then be written as j tt > = j{2t w — (t + 1')) « 
j(t — t w ). Hence, <ti becomes a function of t — t w only, 
and the precise time of switching the flow becomes irrel- 
evant (as it takes place in the steady state after startup). 
This is confirmed by numerical evaluation and also in 
the simulation, by varying 7^. At t = t w , g\ is noth- 
ing but the value of the startup curve; as t — > 00, this 
contribution vanishes since the functions in the integrand 
will generally decay. Furthermore, at f k t w + 7 c /7, a 
maximum occurs as the relevant strain \ jtt'\ ~ 7c when 
the integral is dominated by large t' . This maximum can- 
cels out the corresponding minimum from an, so that the 
stress overshoot is generically expected to vanish. This is 
shown in Fig. [5] for the case ^ w = 0.8. For t w w 7™ ax /7, 
however, the integral determining a\ will have two rel- 
evant contributions, namely t' « t w as well as t' ss 0, 
since the relevant strain at both these times equals 7 C . 
Hence, <J\(t) in this case displays first a maximum and 
then a minimum, as exemplified in Fig. [5] for the case 
7u, = 0.1. The minimum is responsible for <ri(t) not to 
cancel the stress overshoot contained in o"n(i). Obvi- 
ously, as t w — > 0, the contribution <J\(t) vanishes, so that 
perfect elastic recovery is obtained in this limit. While 
the generalized Maxwell model simplifies the interpre- 
tation, the next sections show that broadly equivalent 
curves follow from the present schematic model, which 
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Figure 6. Linear stress responses G e ff = do"/d7||7- 70 |=o.025 
from simulation, fits of the schematic model, and Maxwell 
model as labeled for various waiting strains -y w . The inset 
gives the strain "/ w — 70 required for the stress to reduce to 
zero after shear reversal. To include the simulation point for 
7^ = 75, this value has been divided by 200, using that the 
precise value of -y m is irrelevant in the steady state. 



provides a fundamental basis for this simplification. It 
exhibits a two step process, where the final decay recov- 
ers the properties of Eq. (|T5|) . 

We now turn to a discussion of the effective shear mod- 
uli, obtained as G e g = rfcr/d7|| 7 _ 7o | = o.025- They measure 
the effective elasticity remaining until the shear-driven 
relaxation sets in. Figure [6] displays the results from 
the schematic MCT model (red and blue stars) together 
with the values extracted from the simulation data (blue 
crosses). The agreement is quite satisfactory, noting that 
only the initial value at "f w = is a result of a fitting pro- 
cedure. One observes a marked decrease of G e g around 
7u, ~ 0.1, the position of the overshoot in the startup 
curve. Hence, pre-shear indeed softens the material, but 
only by plastic deformation. Also shown in Fig. [6] arc 
the results from the generalized Maxwell model; quali- 
tatively the same trend is seen. In the Maxwell model, 
we could not find a closed expression to determine 70; 
but the effective shear modulus defined as 7 -> j w for 
t t w + can be evaluated analytically. It again follows 
the trend seen in the figure, albeit with an even lower 
limiting value at large 7^. 

A missing information on the simulated stress-strain 
curve remains 70, the strain value where the stress- free 
state is achieved after strain reversal. Its qualitative be- 
havior can be understood easily. Assuming that G e ff is 
unchanged for small j w , the relations a(t w ) = G e fCf w 
and = <r(t ) = a(t w ) — G c fr(7o — 7™) hold. As a result, 
70 = 2^j w . Linearity holds even if G c g changes, but the 
prefactor is somewhat different, as observed in the simu- 
lations. For large j w , = a(t ) = o- ss -G c fr(7o-7«)) ! and 
the constant value in Fig. [S] is indeed roughly <7 SS /G e g « 
0.4/6 w 0.07. 
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Figure 7. View on the t— t'-plane, with the three different time 
regimes for the density correlator $. As the /i-function de- 
pends only on the absolute value of 7, the one-time-solutions 
$(+^ and $ (_ ^' are equal. 



C. Discussion of the history dependence in the schematic 
model 

To understand the history dependence of the stress- 
strain curves within the schematic ITT-MCT model re- 
quires to analyze the time-dependence of the correlator 
which encodes plastic deformations, and to combine it 
with the one of the elastic coefficient v a , which encodes 
strain-dependent anelasticity. The structure of Eq. (J4j) 
and the time-dependence of the accumulated strain sug- 
gests to split the integration domain into three regimes, 
as shown in Fig. [7] It sketches the corresponding inte- 
gration paths in the t — t'-plane leading to the stresses 
shown in Fig. CD using Eq @. Immediately, the strain 
follows from Eq. ([7]) (remembering that t' < t): 

{7 (*-*') for < H < t < t w , 

7 (2t w - (t + t')) for < f < t w < t, (16) 
-7 (*-*') for t w <t', 

which identifies the three regimes. An integration for 
t',t < t w gives the startup stress-strain curve, the in- 
tegration for t' < t w < t gives the stress contribution 
a 1 in Eq. (|14a|) . and the integration for t w < t',t the 
contribution ajj in Eq. (|14b|) . Recall that the correla- 
tors appearing in the MCT model are transient correla- 
tion functions, i.e., they are formed with the equilibrium 
distribution function always, although the full noncqui- 
librium dynamics appears in their propagator. A corre- 
lation function $(i,t') hence contains only information 
on the flow history in the interval [t',i\. In particular, 
the correlation functions defined for < t' < t < t w are 
given by the startup transient solution for -j(t) = +7 
of the model; the same holds true for those defined for 
t w <t' < t, since there the shear rate is constant and the 
direction of shear cannot enter by symmetry. Because 
the strain in these two time regimes depends on the time 
difference only (see Eq. (HHJ)), the correlator $^(i — t') 



also depends on the time difference only. This greatly 
reduces the complexity of Eq. (J4j) and allows the use of 
well-established integration schemes to solve 



$(7) ( s ) + r $(7) ( s ) + r / da' m W (s - s')^ (s') = 


(17) 

Here, s = t — t' abbreviates the time difference. The 
correlation function for t' < t w < t however remains 
fully two-time dependent, as here the strain depends on 
t + t'. It is convenient to split the integral appearing in 
Eq. (gj) explicitly at t" = t w , 

tw 

t') + r*(t, t') + r [ dt" m(t, t", t')^ {t" - 1') 



+r / dt"m^>(t - t")d t „$(t",t w ) = , 

(18) 

Here, rw^ is the same memory kernel also appearing in 
Eq. ()17|) . Hence, in this formulation, one-time correlators 
or memory kernels previously calculated enter most suc- 
cinctly. Supplied with suitable initial conditions, Eq. (|T8j) 
is a convenient way to deal with instantaneous shear-rate 
switches numerically^. 

Evaluation of a stress (7(7 = f*dt'j(t')) requires in- 
tegrations of the transient shear modulus along lines at 
constant t in the plane shown in Fig. [7] The startup 
curve is obtained along lines of the type from A to B 
before the shear reversal, where all functions depend on 
t—t'. Figures [8] and [9] include the transient density corre- 
lators, the transient shear moduli and the integrands of 
E q. (CD) for this case. As discussed in detail in Refs. [4l|, 
M, l27j|. the linear regime in the stress-strain curve is con- 
nected to the plateau on intermediate times in &(t) and 
G(t), and the stress-overshoot to the negative dip in G(t) 
at late times. Reversing the direction of the shear at a 
late time when the stationary state was reached at 7^, 
correlators and shear moduli change as shown in Fig. [8] 
Three different final states are chosen along the reversing 
stress-strain curve: One in the elastic regime, one in the 
regime where the stress-maximum is now missing, and 
one in the late steady regime, where the steady state is 
reached again. Obviously, the last curves agree with the 
startup curves (except for the trivial minus sign in the 
integrand), as the same transient into a steady state is 
probed. Here the integration path from C to D, where the 
shear rate takes one definite (negative) value, dominates 
and the time region, where the shear rate was different, 
affects the very final transients, only. Moving the final 
time t to correspond to a strain 7 of around 10%, where 
the maximum in the startup curve appears, the transient 
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Figure 8. (a) Shear moduli G(t, t') after shear reversal at 7^ 
in the steady state, shown for three different fixed times t 
(respectively strains 7), as functions of t — t'. Curves marked 
with violet triangle, disc, and diamond symbols correspond 
to points marked with the same symbol in Fig. [1] For times 
t—t' < t w , all curves collapse by construction onto the startup 
result G(t, 0) shown in red (marked with a square). The inset 
shows the corresponding correlators $(t,t') (they overlap on 
the resolution of the figure), (b) Corresponding stress inte- 
grands ±7(t')G(t, t') on a linear t — t' axis. Curves after shear 
reversal which would be negative for t — t' — I are mirrored 
to positive initial decay, where they overlap with the startup 
curve. 



shear modulus is affected in the negative region. There 
the strong variation in v a (t,t') of Eq. <j3j> around 7* be- 
comes important. Along the integration path from D' to 
E', the memory of the +7 shear enters, which causes the 
sign change in G(t, t'); see the curve marked with a circle 
in panel (a) of Fig. [8j The integrand of the generalized 
Green-Kubo relation Eq. (p} varies around zero, see panel 
(b) of Fig. [51 in the region where the stress overshoot 
arose. This leads to a cancellation in the integral, and an 
absence of the maximum in the stress. Shifting the final 
time t to earlier values in the region, where the startup 
stress varies linearly with strain, causes the sign change 
in G(t, t') to move to earlier times t' as well; see the curve 
marked with a triangle in panel (a) of Fig. [5] A stronger 
memory of the shear rate with opposite sign remains. As 
this arises especially during the final relaxation, the re- 
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Figure 9. Like Fig. [8] but now for shear reversal at 7^) in the 
elastic regime. The curves end at a time t — t' = 7^/7, where 
the memory has not decayed to zero yet. This is as expected, 
because the system still responds elastic like before the early 
reversal at 7^; see the blue curve in Fig. [Tj 



gion in G(t, t') where stresses are anticorrelated, the over- 
all integral over the transient shear modulus is smaller. 
As t, respectively 7, lies in the region, where the linear 
increase in the stress occurs, this explains the softening 
of the effective elastic constant G e g; see the linear slope 
indicated in Fig. [T] at the violet triangle. Importantly, in 
the ITT-MCT approach it arises not by a weakening of 
the plateau in the transient correlation functions, rather 
from the aftereffect of the previously stored stresses ac- 
cumulated during the flow with opposite (positive) shear 
rate. It thus originates from the final relaxation and docs 
not hold in the proper limit of infinitesimal strain. We 
chose 7 = 7*/4 to measure G e g. 

The equivalent discussion holds for reversing the flow 
during the elastic regime of the startup stress-strain 
curve. This is indicated by a blue curve in Fig. [T] and 
the corresponding transient functions are shown in Fig. [9] 
Again a triangle, circle, and diamond marks final states 
on the reversing stress-strain curve in the elastic, over- 
shoot, and steady regime. The integration paths in the 
t, t'-plane of Fig. [7] do not change qualitatively, but the 
time t w of switching is now so short that the startup 
correlators along the line A to B have not decayed to 
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Figure 10. In panel (a), linear stress response G c ft = 

da/dyl^ ro i=o.D25 as function of waiting strain -y w in a fluid 

and a glass state; for the latter, two different shear rates are 
shown. The quiescent elastic shear constant Goo is observed 
in startup flow. 

In panel (b), the corresponding curves are given for the rela- 
tive stress overshoots a max /a s — I. 



zero. The system responds clastically, and the correlator 
and the modulus G^' are still of the order of the 
plateau value, / respectively Grx,. The latter value de- 
termines the linear increase of the startup stress-strain 
curve. The complete discussion of the changes in G(t, t') 
done in respect to Fig. [8] carries over to Fig. [9] with the 
sole difference that the two-time region D to E (and D' to 
E') dominates at first. By continuity, it reproduces the 
stress cr(7^|), where the flow reversal took place. Because 
the strain increases with t + t' in this region according to 
Eq. (|16p . this contribution decays rapidly along the re- 
versing stress-strain curve. It may be neglected beyond 
7-7tu, where the initial elastic stress is destroyed. The 
remaining contribution along the lines C to D or C to 
D' is equivalent to the startup flow along A to B, and 
thus after shear rate reversal in the elastic regime, the 
stress-strain curve agrees closely with the initial startup 
curve. 

The preceding discussion of the Bauschinger effect at 
two typical waiting or reversal times t w corresponding to 
strains in the elastic and steady regime enables one to 
rationalize the stress-strain curves for all pre-shears j w . 
Figure [TU] summarizes the main findings for two different 



temperatures close to the glass transition in terms of the 
effective shear constant G e g and the relative amplitude of 
the stress-overshoot <j max ja a — \. Both change from their 
quiescent values as soon as y w becomes of the order of 
7* , where the transient shear modulus captures anticorre- 
latcd stress fluctuations. The variation of G e ff saturates 
when j w > 7**, because then stress correlations have 
become negligible; then also the overshoot o max j a s — 1 
is gone. Increasing the static correlations in the system, 
mimicked by e in the schematic model, increases the elas- 
ticity. Because bare Peclet number Peo -C 1 are con- 
sidered, no qualitative difference holds for a liquid and 
a glassy state. Decreasing the shear rate in the glass, 
the discussed effects remain and just become somewhat 
smaller. In a fluid state, the linear response regime would 
be approached as soon as the dressed Peclet or Wcis- 
senberg number Pc= A fT becomes of order unity. This 
situation where no overshoot arises^ and thus flow rever- 
sal only mildly affects the stress-strain curve is not shown 
here. 



D. Simulation results on average particle motion 

To relate the macroscopic response discussed so far 
to the microscopic dynamics, we turn to a discus- 
sion of simulation results for the density correlation 
functions under shear. Figure [11] shows as exemplary 
cases the self-intermediate scattering functions (tagged- 
particle density-correlation functions) F s (q,t,t') in the 
quiescent equilibrium (t' — > — oo, dashed), under startup 
flow (t' = 0, dash-dotted), and after flow reversal from 
the steady state (f = 7^/7 using 7^ = 0.4, solid lines). 
Each curve has been sampled by averaging over 200 ini- 
tial configurations, and we checked that stationarity is 
achieved comparing with reversal at 7^ = 75. To avoid 
a discussion of shear-advection effects, wave vectors were 
chosen with zero component in the shear direction. The 
equilibrium curves exemplify the typical two-step relax- 
ation found in glass-forming liquids: after an initial fast 
relaxation to an intermediate plateau, the final decay to 
zero is characterized by a large time scale r 10 3 . This 
final decay is much more stretched than exponential re- 
laxation. The plateau value depends on the wave number 
q, and for the correlation functions characterizing tagged- 
particle motion generically decreases with increasing q. 
Steady shear accelerates the dynamics, and induces re- 
laxation from the plateau on a time scale set by I/7 (with 
a <j-dependent prefactor). This relaxation is much closer 
to exponential, and, for the transient correlation function 
shown here, even slightly compressed. These qualitative 
features confirm earlier observation o 6 i 30 i 42 . By symme- 
try, the same is true for the correlation functions prob- 
ing the reversed- flow regime (solid lines in the figure). 
However, as in the macroscopic response, a pronounced 
history dependence is seen: the intermediate plateau is 
much less pronounced, and the final shear-induced decay 
shows much more pronounced stretching, in particular at 
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Figure 11. Simulation results of the self intermediate scat- 
tering functions F s (q,t) for the B-species in the two di- 
rections perpendicular to shear during startup and after 
shear reversal in the stationary state, for wave vectors q = 
[tv/uaa, 2-k/uaa, 3ir/<JAA,4:ir/aAA] in descending order. A to- 
tal 200 independent simulational runs are averaged to obtain 
these graphics. The corresponding equilibrium functions are 
included (dashed lines, label EQ), shifted to coincide at short 
times. 



large q. 

In discussing these dynamics, one has to emphasize the 
difference of these measurable quantities to the transient 
correlation functions used in MCT. The theory defines 
correlation functions as averages over the equilibrium dis- 
tribution; they have the advantage to be better suited 
for subsequent approximations. In the simulation, such 
correlation functions can only be measured (beyond the 
quiescent equilibrium) under startup of steady shear from 
an equilibrated configuration. This is the case shown in 
Fig.Qj]by the dash-dotted lines, which can be compared 
qualitatively — in view of the lack of a g-dependence of 
the latter — to correlators from the schematic model in- 
cluded in the inset of Fig. [8] For the case of flow reversal 
(solid lines), the simulation implies averaging over the 
steady-state distribution corresponding to forward flow, 
and not the quiescent one. We observe that this quan- 
tity shows changes already in the relaxation around the 
plateau, and not just for features of the final relaxation. 

A related quantity that can be intuitively interpreted 
is the mean-squared displacement (MSD), Sr 2 (t,t'), of a 
single particle. We choose randomly a particle of species 
B, and show the results in Fig. [12] Here we show, in ad- 
dition to the cases t' -C 0, t' = 0, and t' = t%, discussed 
above, also the MSDs obtained for t' = and t' = t™ ax . 
For simplicity, we only record the motion in the neutral 
direction; the MSD is then defined as Sr 2 (t) = 3Sz 2 (t) to 
match the spatially averaged MSD of the quiescent sys- 
tem at short times. 8z 2 {t) is expected to become diffusive 
at long times, while the other components of the MSD 
under shear attain supcrdiffusive asymptotes due to (the 
high-density analog of) Taylor dispersion^. 




\i\(t-t') 




Figure 12. (a) Mean-squared displacement of B particles in 
the vorticity direction. Shown are the equilibrium curve for 
T = 0.4 (label EQ), and transient MSDs Sz 2 (t,t') as a func- 
tion of |7|(t — t'), where t' is the time of initial flow start-up 
from equilibrium (EQ) respectively flow reversal according to 
the cases discussed in Fig. [5] (steady-state, s; elastic tran- 
sient, el; strain corresponding to the overshoot, max). The 
equilibrium curve is shifted to coincide at short times, (b) 
Logarithmic derivative /tt t /(t) = rflog 5z 2 (t, t')/d\og t for the 
cases shown in (a). 



It has been established by computer simulation and 
confocal microscopy studies^^ that the stress-overshoot 
in startup flow is accompanied by superdifFusive particle 
motion on intermediate time scales. This is confirmed 
by the startup case shown in Fig. 1121 while the equilib- 
rium MSD remains subdiffusive at all times, the transient 
MSD obtained for t' = follows the equilibrium one up to 
a time that corresponds to a few percent strain. Then, 
it quickly crosses over to the diffusive regime obtained 
under steady-state flow. An effective exponent can be 
assigned to the MSD curves by taking the logarithmic 
derivative, ji(t) = d\og5r 2 (t) / d\ogt. This quantity is 
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shown in the lower panel of Fig. [T2j Since the short-time 
motion of the simulation is ballistic, n(t) = 2 holds for 
t —> 0; in the plateau region, fi(t) drops to values close 
to zero, and for long times, fj,(t) = 1 indicates diffusive 
motion. For times where the fast increase in the if = 
MSD is seen, fi{t) rj 1.3 holds, in qualitative agreement 
with the simulation of Ref . [f| . 

Remarkably, the correlation between stress overshoot 
and superdiffusive motion of the individual particles 
holds beyond the case of startup. In Fig. [T^J we observe 
superdiffusive MSDs for those cases, where the corre- 
sponding stress-strain curve in Fig.[3]shows an overshoot; 
for the case of flow reversal in the steady state, where the 
stress overshoot has vanished, also only subdiffusive mo- 
tion is seen in the MSD. This hints at a strong connection 
between the two phenomena. Such a connection can be 
rationalized^ by invoking a generalized Stokes-Einstein 
relation: the mean-squared displacement is determined 
by a Mori-Zwanzig equation similar to Eq. (j2J). The rel- 
evant memory kernel can be approximated by the gen- 
eralized time-dependent shear modulus appearing in the 
Green-Kubo relation for the stress, Eq. (fTJ). Although it 
does not provide a sufficient condition, one sees that a 
memory kernel with negative portions in the MSD equa- 
tions may lead to superdiffusive behavior. Both stress- 
overshoot and superdiffusive MSD would then arise in 
the same time window characterized by the strain 7*. 



V. CONCLUSIONS 

We have discussed the history-dependent response of 
sheared glass-forming liquids subject to shear flow whose 
direction is instantaneously reversed. Starting from the 
equilibrated liquid, application of a constant shear rate 
7 causes the resulting stress er to increase from zero to 
a steady-state value, undergoing an intermediate maxi- 
mum, the stress overshoot. Flow reversal, or equivalently, 
application of a constant shear rate to a system that is 
in steady state corresponding to the flow in opposite di- 
rection, results in a stress-versus-time curve that has no 
maximum. At the same time, the (effective) elastic shear 
modulus found from the initial rate of change of the stress 
is lowered, and the regime in applied strain shrinks over 
which the behavior is described by this initial linear elas- 
tic regime. These coincident phenomena can be seen as 
analogous to the Bauschinger effect in simple shear flow. 
It should be pointed out that the original Bauschinger 
effect was discussed for metals and steels under com- 
pressional and tensile load. However, given the generality 
of history-dependent flow phenomena in amorphous sys- 
tems, one may expect that the underlying mechanisms 
are related to our discussion. 

We have presented a systematic investigation of the 
gradual change between the two cases of starting shear 
from a quiescent equilibrium and a sheared steady state 
that belongs to the opposite flow. Reversing the flow 
direction after an initial deformation that remains in the 



linear elastic regime, no change in the stress overshoot 
or the elastic coefficients is seen. The position of the 
overshoot indeed, as envisaged already by Bauschinger—, 
marks the end of the reversible elastic regime: at this 
point, plastic deformation gradually takes over, and the 
effective elastic coefficients drop over a small window in 
applied pre-strain. 

Our theoretical analysis follows the integration- 
through-transicnts (ITT) scheme upon which the mode- 
coupling theory for colloidal rheology is based. The cen- 
tral physical mechanism expressed in our model is one of 
temporal history dependence: glass-forming liquids are 
visco-elastic, with a large time window r over which the 
fate of past density fluctuations affects the present re- 
sponse. It is this history dependence that causes the 
Bauschinger effect and similar pre-shear-dependent phe- 
nomena in the model, rather than spatial anisotropics 
induced by the pre-shear (since no information on spa- 
tial variation of fluctuations is kept). This is in align- 
ment with the microscopic MCT and molecular-dynamics 
simulations of the sheared thermal flui d 30 i 43 , where it 
was found that the steady state exhibits almost isotropic 
(one-time) pair correlation functions. This may be differ- 
ent in the athermal limit. It does, however, qualitatively 
explain the disappearance of the stress overshoot, and 
our findings for the tagged-particle correlation functions 
and mean-squared displacements. Immediately, this dis- 
cussion predicts, that inserting a waiting time between 
the opposing flows would cause the Bauschinger effect 
to weaken and go away. Any quiescent waiting period 
before shear-reversal would enable relaxation to proceed 
and would ultimately lead to startup curves from the qui- 
escent state. When shear is started from the quiescent 
equilibrium configuration, nearest-neighbor cages need to 
be broken before the accumulated strain becomes effec- 
tive in enhancing the relaxation dynamics. The critical 
strain is hence connected to a typical cage size; some 
percent of the particle diameter following a Lindemann- 
type criterion for melting. The sheared steady state is 
then characterized by less strong cages, and the direc- 
tionality of the flow does not play a major role in this 
weakening effect. Hence, subsequent shear in the oppo- 
site direction will lead to an earlier and more gradual re- 
laxation of correlation functions (cf. Fig. [TTj) . The mean- 
squared displacements, observed to follow the unsheared 
equilibrium curve up to times corresponding to the criti- 
cal strain in the case of startup shear, remain much closer 
to the sheared steady-state ones in the case of flow re- 
versal from such a steady state. As a consequence, the 
sharp upturn observed for the startup case, connected to 
superdiffusive motion and a stress overshoot, vanishes. 

At present, the analysis of the initial stress overshoot 
after startup of steady shear is based on a schematic 
model that does have a number of fit parameters whose 
numerical ratio has to be determined phenomenologi- 
cally. However, the model is rooted in an understanding 
of the microscopic MCT—. In simplifying the latter to ar- 
rive at a schematic model, one aims to keep those essen- 



14 



tial features of the solutions qualitatively alike that the 
model shall describe; in this sense, the model used here is 
a minimal model containing stress-overshoot phenomena. 
It is important to note that the subsequent description 
of the fate of the stress overshoot under time- varying 
flow drops out naturally from the equations of motion. 
In this sense, the comparisons presented in Figs. and 
[3] are qualitatively different - the first can be viewed 
as a model-motivated paramctrization of the simulation 
data. The second - the comparison with the data on the 
Bauschinger effect - entails a parameter-free theoretical 
prediction. 
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